Dynamics of collapsing and exploding Bose-Einstein condensate 



Sadhan K. Adhikari 

Instituto de Fisica Teorica, Universidade Estadual Paulista, 01.405-900 Sao Paulo, Sao Paulo, Brazil 

(February 1, 2008) 



Recently, Donley et al. performed an experiment on the dynamics of collapsing and exploding 
Bose-Einstein condensates by suddenly changing the scattering length of atomic interaction to a 
large negative value on a preformed repulsive condensate of 85 Rb atoms in an axially symmetric 
trap. Consequently, the condensate collapses and ejects atoms via explosions. We show that the 
accurate numerical solution of the time-dependent Gross-Pitaevskii equation with axial symmetry 
can explain some aspects of the dynamics of the collapsing condensate. 



Since the successful detection |lj |S|] of Bose-Einstein 
condensates (BEC) in dilute bosonic atoms employing 
magnetic trap at ultra-low temperature, one problem of 
extreme interest is the dynamical study of the formation 
and decay of BEC for attractive atomic interaction . 

For attractive interaction the condensate is stable for 
a maximum critical number iVcr of atoms 0] . When the 
number of atoms increases beyond this critical number, 
due to interatomic attraction the radius of BEC tends 
to zero and the central density tends to infinity. Con- 
sequently, the condensate collapses emitting atoms until 
the number of atoms is reduced below iVcr and a sta- 
ble configuration is reached. With a supply of atoms 
from an external source the condensate can grow again 
and thus a series of collapses can take place, which was 
observed experimentally in the BEC of 7 Li with attrac- 
tive interaction ||. Theoretical analyses based on the 
time-dependent mean-field Gross-Pitaevskii (GP) equa- 
tion also confirm the collapse [p[-pl^. 

It is possible to manipulate the interatomic interaction 
by an external magnetic field via a Feshbach resonance 
@. Roberts et al. @ and Cornish ct al. @] (at JILA) 
exploited this idea to suddenly change the atomic scat- 
tering length by a large amount in experiment. They 
have been able to even change the sign of the scatter- 
ing length, thus changing a repulsive condensate in to 
an attractive one and vice versa. Consequently, a sta- 
ble preformed repulsive condensate is suddenly turned 
to a highly explosive and collapsing attractive conden- 
sate. In a classic experiment performed at JILA Donley 
et al. jL5j studied the dynamics of collapsing and ex- 
ploding condensates formed of 85 Rb atoms. The natural 
scattering length of 85 Rb atoms is negative (attractive). 
By exploiting the Feshbach resonance they made it pos- 
itive (repulsive) in the initial state. They have provided 
a quantitative measurement of this explosion by count- 
ing the number of emitted and remaining atoms in the 
condensate as a function of time until an equilibrium is 
reached. They also measured the energy distribution of 
the emitted atoms. They claim that their experiment 
reveal many interesting phenomena that challenge theo- 
retical models. 

We demonstrate that some aspects of the above col- 
lapse and explosion of the attractive condensate of 85 Rb 



atoms can be explained from an accurate numerical solu- 
tion of the GP equation in an axially symmetric trap, 
where we include a quintic three-body nonlinear re- 
combination term which accounts for the decay of the 
strongly attractive condensate. The numerical method, 
we use, for the solution of the time-dependent GP equa- 
tion with an axially symmetric trap has appeared else- 
where @j22|. 

There have been other theoretical studies |]rj|-|l0|| to deal 
with dynamical collapse including an absorptive term to 
account for the loss of particles. Instead of attempting 
a full numerical solution of the GP equation with axial 
symmetry, these investigations used various approxima- 
tions to study the time evolution of the condensate or 
employed a spherically symmetric trap. 

The time-dependent Bose-Einstein condensate wave 
function ^(r; r) at position r and time r may be de- 
scribed by the following mean-field nonlinear GP equa- 
tion HI 
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Here m is the mass and N the number of atoms in the 
condensate, g = ATiTi 2 a/m the strength of interatomic in- 
teraction, with a the atomic scattering length. The trap 
potential with cylindrical symmetry may be written as 
V{v) = ^mui 2 (r 2 +X 2 z 2 ) where u> is the angular frequency 
in the radial direction r and Xlo that in the axial direc- 
tion z. We are using the cylindrical coordinate system 
r = (r, 8, z) with 9 the azimuthal angle. The normaliza- 
tion condition of the wave function is / dr|\l/(r; r)| 2 = 1. 

In the absence of angular momentum the wave function 
has the form ^r; r) = ip(r, z; r). Now transforming to di- 
mensionless variables defined by x — \jlrjl, y — \/2z/l, 
t = toj, I = yh] (mw), and 
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where n = Na/l. The normalization condition of the 
wave function becomes 
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A/norm = 2% / dx dy\ip(x,y;t)\ 2 x = 1. (4) 
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The root mean square (rms) sizes Xrms and Jfrms are 
defined by 

Zrms = -A/norm 271 " / ^ / dy\tp(x,y;t)\ 2 x, (5) 
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yrms = Miorm 27r / da; / dy|</j(:E, y; t)| V* X - (6) 



It is now appropriate to calculate the parameters of the 
present dimensionless GP equation (^) corresponding to 
the experiment at JIL A fl5|| . We follow the notation and 
nomenclature of Ref. [15fl. Their radial and axial trap 
frequencies are ^ racuai i = 17.5 Hz and f ax i a i = 6.8 Hz, 
respectively, leading to A = 0.389. The harmonic oscil- 
lator length I for 85 Rb atoms for u> — 2tt x 17.5 Hz is 
/ = y^K/Jmco) = 25905 A. One unit of time t of Eq. 
(|) is 1/w or 0.009095 s. They prepared a stable 85 Rb 
condensate of No = 16000 atoms with scattering length 
a = ainitial = 7ao, ao = 0.5292 A, such that the initial 
n = 2.288. Then during an interval of time 0.1 ms the 
scattering length was ramped to a = a co n a p Se = — 30ao 
such that final n = —9.805. The final condensate is 
strongly attractive and unstable and undergoes a se- 
quence of collapse and explosion. 

The sequence of collapse in many theoretical studies 
has been explained by introducing an absorptive imag- 
inary three-body quintic interaction term of strength £ 
responsible for recombination loss. Consequently Eq. (j^) 
becomes 01 
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There are many ways to account for the loss mechanism 
HQ- It is quite impossible to include them all in a self 
consistent fashion. Here we simulate the loss via the most 
important quintic three-body term with parameter £. We 
use this single parameter to reproduce the experiment at 
JILA |J). For £^ in Eq. (fi), AWn + 1. 

We solve the GP equation (7) numerically by time iter- 
ation with a given initial solution. In the time-evolution 
of the GP equation the radial and axial variables are 
dealt with in independent steps. For this purpose we 
discretize it using time step A = 0.001 and space step 
0.1 for both x and y spanning x from to 15 and y 
from —25 to 25. This domain of space was sufficient to 



encompass the whole condensate wave function even dur- 
ing and after explosion and collapse. The preparation of 
the initial repulsive wave function is now a routine job 
and was done by increasing the nonlinearity n of the GP 
equation (Q) by 0.0001 in each time step A during time 
iteration starting with the known harmonic oscillator so- 
lution of Eq. (^) for n = £ = | flq |. The initial value 
of n(= 2.288) was attained after 22880 time steps. The 
nonlinearity n is then ramped from 2.288 to —9.805 in 
0.1 ms. As one unit of dimensionless time t is 0.009095 s, 
0.1 ms corresponds to 11 steps of time A. In the present 
simulation, n is ramped from 2.288 to —9.805 in the GP 
equation by equal amount in 11 steps. The absorptive 
term £ was set equal to zero during above time iteration. 
Now the system is prepared for the simulation of collapse 
and explosion. 

For the remaining simulation the nonlinear term is 
maintained constant and a nonzero value of £ is chosen. 
The time-evolution of the GP equation is continued as 
a function of time t = T evo i ve starting at 0. Because 
of the nonzero dissipative term £, the normalization (^) 
is no longer maintained, and the number of remaining 
atoms A^ in the condensate is given by A^ = AoA/norm, 
where Nq is the initial number. The time-evolution is 
continued using time step A = 0.001. After a small ex- 
perimentation it is found that £ = 3.12 fits the data of 
the experiment at JILA for a C0 H a pse = ~ 30ao — their 
Fig. 1(b). This value of £ was used in all simulations 
reported here. The remaining number of atoms vs. time 
is plotted in Fig. 1. 

It is pertinent to compare this value of £(= 3.12) with 
the experimental measurement of three-body loss rate on 
85 Rb |p3[ . For that, we need to restore the proper dimen- 
sions in the three-body term of Eq. (^) by rewriting it 
as i£,(aN/l) 2 \ip\ 4 l 6 /8 and equating it to the convensional 
form iK 3 N 2 \?P\ 4 /2Lu, where K 3 = (4.24t°;™ ± 0.85) x 
10~ 25 cm 6 /s is the experimental rate [p3[ . From this 
we find K 3 = £a 2 Z 4 a;/4. Under experimental condition 
of an external magnetic field of 250 gauss on 85 Rb |2l| 
the scattering length was a ~ — 350ao. Consequently, the 
present value of £(= 3.12) corresponds to K 3 ~ 13xl0~ 25 
cm 6 /s, about three times larger than the experimental 
(7) value above. 

In the experiment at JILA [jl5| it was observed that the 
strongly attractive condensate after preparation remains 
stable with a constant number of atoms for an interval of 
time ^collapse' This behavior is physically expected. Im- 
mediately after the jump in scattering length from 7ao to 
— 30ao, the attractive condensate shrinks in size during 
''collapse un ^ u t ne central density increases to a maxi- 
mum. Then the absorptive three-body term takes full 
control to initiate the sequence of collapse and explosion. 
Consequently, the number of atoms remains constant for 
T evolve < 'collapse- ^ ne P resen t result (full line) also 
shows a similar behavior. However, in this simulation 
the absorptive term is operative from T evo i ve = and 
the atom number decreases right from beginning, albeit 
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at a much smaller rate for 7" evo i ve < Collapse' 

Donley et al. |0| repeated their experiment with differ- 
ent values of final scattering length. Using the same value 
of £(= 3.12), we also repeated our calculation with the 
final scattering lengths: G caua p Se = — 6.7ao and — 250ao- 
These results are also plotted in Fig. 1. The initial de- 
lay i co n a p Se in starting the explosion is large for small 
l a collapsel as we see * n ^ After the initial delay, the 
decay rates of the number of atoms are quite similar as 
observed in the experiment. Similar effect was observed 
in the experiment for an initial condensate of 6000 atoms 
as shown in their Fig. 2. 

After a sequence of collapse and explosion, Donley et 
al. [ fj"5"[ observed a "remnant" condensate of ^remnant 
atoms at large times containing a certain fraction of ini- 
tial No atoms. Experimentally, the fraction of atoms 
that went into the remnant decreased with | a collapse I 
and was ~ 40% for |a co n a p ge | < 10a and was ~ 10% 
l a collapsel > lOOao Figure 1 also shows this behavior. 
We repeated the simulation of Fig. 1 with A*o = 6000 
where we also found a similar behavior. The number in 
the remnant condensate could be much larger than TVcr- 
The above evolution of the condensate after the jump 
in scattering length to — 30ao for No — 16000 can be 
understood from a study of the wave function and we 
display the central part of the wave function in Fig. 2 
for T evo i ve = 0,3.5,4, and 10 ms. The wave function 
immediately after jump at time T evo j ve = is essentially 
the same as that before the jump at —0.1 ms. There is 
not enough time for the wave function to get modified 
at T evo j ve = 0. From Fig. 2 we find that at 3.5 ms the 
wave function is only slighty narrower than at ms but 
still smooth and has not yet collapsed significantly. As 
T evolve increases, the wave function contracts further and 
the explosion starts. At 4 ms some spikes have appeared 
in the wave function showing the beginning of explosions 
and loss. From the study of the wavefunctions we find 
that the explosions start at T evo j ve = £ co llapse — 3-7 nis 
in agreement with the experiment at JILA. We also find 
that at 3.7 ms before the loss began the bulk BEC did 
not contract dramatically as also observed in the experi- 
ment. In numerical simulation for this case we find that 
at T evo i ve = 0,a;rms = 2.98/im and yrms = 4.21/mi and 
at r evo j ve = 3.7 ns, x rms = 2.53/Ltm and y r ms = 4.10/zm. 
From Fig. 2 we see that at 10 ms the wave function is 
very spiky corresponding to many violent ongoing explo- 
sions. 

Donley et al. |l5| also observed that the remnant con- 
densate oscillated in a highly excited collective state 
with approximate frequencies 2^ ax j a ^ and 2f ra( jj a j be- 
ing predominantly excited. The measured frequencies 
were 13.6(6) Hz and 33.4(3) Hz. To find if this behavior 
emerges from the present simulation we plot in Fig. 3 
sizes Xrms and yrms vs. time for the condensate after 
the jump in the scattering length to — 30ao from 7ao for 
Nq = 16000. Excluding the first 20 ms when the remnant 



condensate is being formed, we find a periodic oscillation 
in xrms and yrms with frequencies 34 Hz and 13.6 Hz, 
respectively. In addition, for xrmsj superposed on 34 
Hz we note a prominent frequency of about 3.5 Hz, e.g., 

°- 2l/ radial- 

Though we have explained some aspects of the exper- 
iment at JILA, there are certain detailed features that 
have not been obtained in this Letter. Donley et al. 
Jl5[ have classified the emitted atoms in three categories: 
burst, jet and missing atoms. Detailed study of the wave 
function is needed to identify the jet atoms that appear 
possibly from the spikes in the wave function when the 
collapse is suddenly interrupted. The jets were also not 
always found to be symmetric about the condensate axis 
|]l5f . This makes an axially symmetric mean-field model 
as used in this Letter inappropriate for a description of 
the jet atoms. Moreover a clear-cut distinction between 
the burst and missing atoms emitted during explosions 
seems to be difficult in the present model. Also because 
of the missing (undetected) atoms it is difficult to predict 
the energy distribution of the burst atoms in a mean-field 
analysis which will only yield the total energy distribu- 
tion of the burst plus missing atoms. A careful analysis 
of the energy of the emitted atoms is required for ex- 
plaining these detailed features and would be a welcome 
future theoretical investigation. 

In conclusion, we have employed a numerical simula- 
tion based on the accurate solution Jl(| of the mean-field 
Gross-Pitaevskii equation with a cylindrical trap to study 
the dynamics of collapse and explosion as observed in the 
recent experiment at JILA In the GP equation we 
include a quintic three-body nonlinear recombination loss 
term which accounts for the decay of the strongly attrac- 
tive condensate. The results of the present simulation 
accounts for the essential features of that xperiment. We 
find that the exact inclusion of the asymmetry parameter 
A in the simulation is essential for a proper description 
of the experiment. In the experiment at JILA a strongly 
attractive 85 Rb condensate was prepared by ramping the 
scattering length to a large negative value and the sub- 
sequent decay of the condensate was measured. We have 
been able to understand the following features of this dy- 
namics from the present numerical simulation: (1) The 
condensate undergoes a sequence of collapse and explo- 
sion and finally stabilizes to a remnant condensate at 
large times containing about ~ 10% — 40% of the initial 
number No- The number in the remnant condensate can 
be larger than the critical number for collapse Nqt for the 
same atomic interaction. (2) Both in the experiment and 
our simulation the remnant condensate executes radial 
and axial oscillations in a highly excited collective state 
for a long time with frequencies 2f ra( jj a j and ^ ax [ a \- (3) 
The rate of decay of the number of particles is reasonably 
constant independent of Nq and « co n a p Se - (4) After the 
sudden change in the scattering length to a large negative 
value, the condensate needs an interval of time i co vj a p Se 
before it experiences loss via explosions. Consequently, 
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the decay starts after the interval of time i co n a p Se . 
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in ms. The solid circles fll5[ are from experiment with 
a collapse = — 30ao and the theoretical curves are labeled 
by the values of a collapse . 

2. The central part of the dimensionless wave func- 
tion \(f>(x,y)\ = \ip(x,y)/x\ of the condensate on 0.1 x 0.1 
grid after the jump in the scattering length of a BEC 
of 160 85 Rb atoms from 7a to a co n a p Se = — 30ao at 
times T evo j ve = 0, 3.5 ms, 4 ms, and 10 ms. 

3. The dimensionless rms sizes Xrms line) and 
2/rms (dashed line) expressed in units of lj\[2 after the 
jump in the scattering length of a BEC of 160 85 Rb 
atoms from a m iy a l = 7ao to » co u a p SC = — 30ao as func- 
tions of time T evo i ve - 



Figure Caption: 

1. Number of remaining atoms in the 85 Rb conden- 
sate of 16000 atoms after ramping the scattering length 
from a m jti a l = 7a to a co n a p Se = — 6.7a ,— 30a and 
-250a in 0.1 ms as a function of evolution time 7"gyolve 
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